TABLEAU EXERCISE

Mapping of toto outlets in tableau by using wgs coordinate format and indicating the X and Y to longtitude and latitude data formats respectively. For project coordinates systems like svy, need to do conversion $Map here$
Realis data does not have X an Y coordinates, but Tableau is able to consume the API from www.onemap.gov.sg to plot out as a map, and then return the latitude and longtitude for you. Realis data is from URA residental transaction website.

Use the smu library CEIC
GADM - get administrative boundary map for free in KMZ or Shapeflie format
data.gov.sg - get boundaries map
singstat - get residents dat by planning subzone

R EXERCISE

Ex 1:Mapping Geospatial Point Data with R

To load packages required today:

packages = c('tidyverse', 'sf', 'tmap')

for(p in packages){
  if(!require(p, character.only = T)){
    install.packages(p)
  }
  library(p, character.only = T)
}

sgpools <- read_csv ("data/aspatial/SGPools_svy21.csv")

head (sgpools, 3)

# A tibble: 3 x 7
  NAME   ADDRESS POSTCODE XCOORD YCOORD `OUTLET TYPE` `Gp1Gp2 Winning~
  <chr>  <chr>      <dbl>  <dbl>  <dbl> <chr>                    <dbl>
1 Livew~ 2 Bayf~    18972 30842. 29599. Branch                       5
2 Livew~ 26 Sen~    98138 26704. 26526. Branch                      11
3 Sport~ Lotus ~   738078 20118. 44888. Branch                       0

combining the xcoord and ycoord into a new column of spatial data and the order of the coordinates must be longtitude-latitude.

sgpools_sf <- st_as_sf(
  sgpools,
  coords = c("XCOORD", "YCOORD"),
          crs=3414)

head (sgpools_sf, 3)

Simple feature collection with 3 features and 5 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 20117.93 ymin: 26525.7 xmax: 30841.56 ymax: 44888.06
Projected CRS: SVY21 / Singapore TM
# A tibble: 3 x 6
  NAME       ADDRESS           POSTCODE `OUTLET TYPE` `Gp1Gp2 Winning~
  <chr>      <chr>                <dbl> <chr>                    <dbl>
1 Livewire ~ 2 Bayfront Avenu~    18972 Branch                       5
2 Livewire ~ 26 Sentosa Gatew~    98138 Branch                      11
3 SportsBuz~ Lotus Lounge, Le~   738078 Branch                       0
# ... with 1 more variable: geometry <POINT [m]>

Using tmap to wrap everything into a leaflet format. tmap works like ggplot2, for map creation. all maps must start with tm_shape

tmap_mode("view")

tm_shape(sgpools_sf) +
tm_bubbles(col = "green",
            size = 0.5,
            border.col = "black",
            border.lwd = 1,
           interactive = TRUE)

Faceting 2 maps for synchronisation view

sync = True will ensure the 2 maps move in tandem

tm_shape(sgpools_sf) +
tm_bubbles(col = "OUTLET TYPE",
            size = "Gp1Gp2 Winnings",
            border.col = "black",
            border.lwd = 0.5,
           interactive = TRUE) +
  tm_facets(by= "OUTLET TYPE",
            nrow = 1,
            sync = TRUE)

Ex 2: Choropleth Mapping with R

use st_read function to convert the .shp file into 1 single tibble object, to get an administrative boundary map, and a column of geometry (string of coordinates)

mpsz <- st_read(dsn = "data/geospatial",
                layer = "MP14_SUBZONE_WEB_PL")

Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `C:\farahfoo\VA\In_Class_ex\In-class_Ex07\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21

head (mpsz, 1)

Simple feature collection with 1 feature and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 30794.28 ymin: 28369.47 xmax: 32362.39 ymax: 30140.01
Projected CRS: SVY21
  OBJECTID SUBZONE_NO    SUBZONE_N SUBZONE_C CA_IND   PLN_AREA_N
1        1          1 MARINA SOUTH    MSSZ01      Y MARINA SOUTH
  PLN_AREA_C       REGION_N REGION_C          INC_CRC FMEL_UPD_D
1         MS CENTRAL REGION       CR 5ED7EB253F99252E 2014-12-05
    X_ADDR   Y_ADDR SHAPE_Leng SHAPE_Area
1 31595.84 29220.19   5267.381    1630379
                        geometry
1 MULTIPOLYGON (((31495.56 30...

Read in the population data

popagsex <- read_csv("data/aspatial/respopagsex2000to2018.csv")

head (popagsex, 2)

# A tibble: 2 x 6
  PA         SZ        AG     Sex       Pop  Time
  <chr>      <chr>     <chr>  <chr>   <dbl> <dbl>
1 Ang Mo Kio Cheng San 0_to_4 Males     810  2000
2 Ang Mo Kio Cheng San 0_to_4 Females   700  2000

prepare the population data for use

spread function is an old function already. Try replacing with pivot_wider
note the difference in the code to sum the Young and Economy active. Young lists all the column names, but economy active counts the column number.

popagsex2018_male <- popagsex %>%
filter(Sex == "Males") %>%
filter(Time == 2018) %>%
spread(AG, Pop) %>%
mutate(YOUNG = `0_to_4`+`5_to_9`+`10_to_14`+
`15_to_19`+`20_to_24`) %>%
mutate(`ECONOMY ACTIVE` = rowSums(.[9:13]) + 
                          rowSums(.[15:17]))%>%
mutate(`AGED`=rowSums(.[18:22])) %>%
mutate(`TOTAL`=rowSums(.[5:22])) %>%
mutate(`DEPENDENCY` = (`YOUNG` + `AGED`)/`ECONOMY ACTIVE`) %>%

# this converts the PA and SZ to upper case to left_join later on
mutate_at(.vars = vars(PA, SZ),.funs = funs(toupper)) %>%
  
select(`PA`, `SZ`, `YOUNG`,`ECONOMY ACTIVE`, `AGED`,`TOTAL`, `DEPENDENCY`) %>%
filter(`ECONOMY ACTIVE` >0)

head (popagsex2018_male, 2)

# A tibble: 2 x 7
  PA         SZ          YOUNG `ECONOMY ACTIVE`  AGED TOTAL DEPENDENCY
  <chr>      <chr>       <dbl>            <dbl> <dbl> <dbl>      <dbl>
1 ANG MO KIO ANG MO KIO~   660             1290   310  2260      0.752
2 ANG MO KIO CHENG SAN    3110             8020  2440 13570      0.692

left join the 2 tables using SUBZONE_N and SZ

glimpse (popagsex2018_male)

Rows: 233
Columns: 7
$ PA               <chr> "ANG MO KIO", "ANG MO KIO", "ANG MO KIO", "~
$ SZ               <chr> "ANG MO KIO TOWN CENTRE", "CHENG SAN", "CHO~
$ YOUNG            <dbl> 660, 3110, 2900, 2370, 890, 1860, 1050, 235~
$ `ECONOMY ACTIVE` <dbl> 1290, 8020, 7610, 6370, 1730, 4800, 2200, 5~
$ AGED             <dbl> 310, 2440, 2710, 2130, 540, 1430, 640, 2020~
$ TOTAL            <dbl> 2260, 13570, 13220, 10870, 3160, 8090, 3890~
$ DEPENDENCY       <dbl> 0.7519380, 0.6920200, 0.7371879, 0.7064364,~

mpsz_agemale2018 <- left_join(mpsz,popagsex2018_male,by = c("SUBZONE_N"="SZ"))

head (mpsz_agemale2018, 2)

Simple feature collection with 2 features and 21 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 28160.23 ymin: 28369.47 xmax: 32362.39 ymax: 30177.73
Projected CRS: SVY21
  OBJECTID SUBZONE_NO    SUBZONE_N SUBZONE_C CA_IND   PLN_AREA_N
1        1          1 MARINA SOUTH    MSSZ01      Y MARINA SOUTH
2        2          1 PEARL'S HILL    OTSZ01      Y       OUTRAM
  PLN_AREA_C       REGION_N REGION_C          INC_CRC FMEL_UPD_D
1         MS CENTRAL REGION       CR 5ED7EB253F99252E 2014-12-05
2         OT CENTRAL REGION       CR 8C7149B9EB32EEFC 2014-12-05
    X_ADDR   Y_ADDR SHAPE_Leng SHAPE_Area     PA YOUNG ECONOMY ACTIVE
1 31595.84 29220.19   5267.381  1630379.3   <NA>    NA             NA
2 28679.06 29782.05   3506.107   559816.2 OUTRAM   620           2210
  AGED TOTAL DEPENDENCY                       geometry
1   NA    NA         NA MULTIPOLYGON (((31495.56 30...
2 1160  3990  0.8054299 MULTIPOLYGON (((29092.28 30...

note that this map has missing values as there are some ares with no population

tmap_mode("plot")

qtm(mpsz_agemale2018,fill ="DEPENDENCY")

Drawing the map with tm_fill and tm_border

tm_shape(mpsz_agemale2018) +
tm_polygons()

tm_shape(mpsz_agemale2018) +
tm_polygons("DEPENDENCY")

lwd = border line width. The default is 1,
alpha= transparency number between 0 (totallytransparent) and 1 (not transparent). By default,the alpha value of the col is used (normally 1),
col= border colour, and
lty= border line type. The default is “solid”.

tm_shape(mpsz_agemale2018) +
tm_fill("DEPENDENCY")

tm_shape(mpsz_agemale2018) +
tm_fill("DEPENDENCY")+
  tm_borders(lwd = 0.1,
             alpha = 1)

Data classification methods of tmap

Most choropleth maps employ some method of dataclassification. The point of classification is to take alarge number of observations and group them into data ranges or classes.

Compare the 2 styles below: quantile vs equal, the intervals are different

tm_shape(mpsz_agemale2018) + 
  tm_fill("DEPENDENCY",n =8,style ="quantile") +
tm_borders(alpha =0.5) +
  tm_layout(main.title ="Quantile interval distribution")

tm_shape(mpsz_agemale2018) + 
  tm_fill("DEPENDENCY",n =8,style ="equal") +
tm_borders(alpha =0.5) +
  tm_layout(main.title ="Equal interval distribution")

adding colour palette

tm_shape(mpsz_agemale2018) +
tm_fill("DEPENDENCY", n = 6,style ="quantile",palette ="Blues") +
  tm_borders(alpha =0.5)

just add a “-” sign infront to inverse the colour scale

tm_shape(mpsz_agemale2018) +
tm_fill("DEPENDENCY", n = 6,style ="quantile",palette ="-Blues") +
  tm_borders(alpha =0.5) +
  tm_layout(main.title ="Reversing colour scheme",
  main.title.position ="center",
main.title.size =1,
legend.height =0.45,
legend.width =0.35,
legend.outside =FALSE,
legend.position = c("right", "bottom"),
frame =FALSE)

Add ing Cartographic Furniture into the tmap

Add in compass, scalebar, gridlines into the map

tm_shape(mpsz_agemale2018)+
tm_fill("DEPENDENCY",style ="quantile",palette ="Blues",title ="No. of persons") +
tm_layout(main.title ="Distribution of Dependency Ratio \nby planning subzone",
main.title.position ="center",main.title.size =1.2,
legend.height =0.45,
legend.width =0.35,
frame =TRUE) +
tm_borders(alpha =0.5) +
tm_compass(type="8star", size =2) +
tm_scale_bar(width =0.15) +
tm_grid(lwd =0.1, alpha =0.2) +
tm_credits("Source: Planning Sub-zone boundary from Urban Redevelopment Authorithy (URA)\n and Population data from Department of Statistics DOS",
             position = c("left","bottom"))

Faceting the maps and custom view to show applicable area only

tm_shape(mpsz_agemale2018) +
  tm_fill("DEPENDENCY",
          style = "quantile",
          palette = "Blues",
          thres.poly = 0) + 
  tm_facets(by="REGION_N",
            free.coords=TRUE,
            drop.shapes=TRUE) +
  tm_layout(legend.show = FALSE,
            title.position = c("center",
                               "center"), 
            title.size = 20) +
  tm_borders(alpha = 0.5)

In class exericse 07